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We propose a two-step protocol for inverting ultrafast nonlinear spectroscopy experiments on a molecular aggregate to 
determine the excited state density matrix at times following laser excitation. The first step is a deconvolution of the 
experimental signal to determine a pump-dependent response function. The second step inverts the quantum state of 
the system from this response function, given a model for how the system evolves following the probe interaction. We 
demonstrate this inversion analytically and numerically for a dimer model system, and evaluate the feasibility of scaling 
it to larger molecular aggregates such as photosynthetic protein-pigment complexes. 



I. INTRODUCTION 

Ultrafast nonlinear spectroscopy allows us to experimen- 
tally observe excited state dynamics in molecular aggregates, 
and in particular, energy transfer essential to the function of 
natural light harvesting systems 1-3 . The existence of these ex- 
perimental tools prompts a natural question: is it possible to 
use spectroscopic measurements to directly infer the excited 
state of such systems? A complete answer to this question 
would be a procedure for quantum state tomography (QST), 
that is, for reconstruction of the full density matrix describ- 
ing the quantum state 4 ' 5 . State tomography is a technique that 
has found widespread application for validating and character- 
izing quantum devices designed as components for quantum 
computation. Such full characterization of an exciton state 
over multiple pigments, beyond a mere classical probability 
distribution, would offer information essential to understand- 
ing the explicitly quantum features of energy transport, which 
include coherence 6 , entanglement 7 and possibly other types 
of non-trivial quantum dynamics 8-11 . In this work, we show 
that under appropriate conditions and assumptions, QST of 
excited states can be performed from the results of a series of 
pump-probe type ultrafast spectroscopies. 

The most sophisticated non-linear technique for resolving 
energy transfer dynamics is the two-dimensional (2D) photon- 
echo, in which the time delays between three ultrafast pulses 
are manipulated to provide a 2D map between pump and probe 
frequencies at fixed time delays 12,13 . These two-dimensional 
maps allow for direct visualization of the relationship between 
excitation and emission energies at a fixed delay time. More 
formally, 2D spectroscopy is usually interpreted in the limit 
of impulsive interactions. In this approximation, it provides 
snapshots of the 3rd-order response function 1 . Important ap- 
plications of 2D spectroscopy in photosynthetic systems have 
included resolving energy transfer pathways 14 and the dynam- 
ics of electronic quantum beats 6,15-17 . In contrast, pump-probe 
spectroscopy (also known as transient absorption) is a simpler 
type of 3rd-order spectroscopy that historically predates 2D. In 
a pump-probe setup, a pump pulse excites the system, which is 
probed at some time later by a probe pulse. Because of its rel- 
ative ease of experimental implementation, pump-probe was 



used to follow ultrafast energy transfer dynamics in photosyn- 
thesis long before 2D spectroscopy. For example, it provided 
the first evidence of electronic quantum beats in photosynthetic 
pigment-protein complexes, in 1997 18 . Pump-probe provides 
less information than the 2D photon echo, because the pump- 
probe signal can be obtained by integrating over the excitation 
axis in a 2D spectra 19 . However, for the purposes of this work, 
pump-probe has a clear advantage, namely, that it can be di- 
rectly interpreted as a measurement of the state created by the 
pump pulse. Formally, the pump dependence is entirely con- 
tained within the change in the density matrix of the system 
after interacting with the pump 20 . 

In the past, time-resolved spectra such as pump-probe have 
been analyzed by simultaneously or concurrently fitting spec- 
tral components, known as decay associated or species associ- 
ated spectra, with a kinetic model 21,22 . When based strictly on 
a series of time-resolved spectra, this procedure is known as 
global analysis; when placed in the context of a specific physi- 
cal model, it is known as target analysis 21 . These techniques 
are powerful, as evidenced by their widespread application to 
experiments. However, much of the kinetic information they 
reveal can be seen more directly in 2D spectra. Moreover, ki- 
netic models, although adequate for many purposes, cannot 
describe more complex dynamics, such as those deriving from 
quantum beats or from a non-Markovian bath. Our approach 
to QST side-steps the issues of extending such integrated anal- 
yses by focusing on identifying the quantum state directly. 

Recently, it was shown that a combination of photon-echo 
measurements of excitonic systems can be combined to per- 
form quantum process tomography of excitonic dimers, either 
by using differently colored pulses 23 or by combining peak 
amplitudes from a set of 2D spectra 24 . Process tomography 4,25 
is more general than state tomography, because it specifies 
the full set of possible quantum evolutions for a system given 
any initial condition. This makes it well suited to character- 
izing gates for quantum computation, which are by definition 
designed to handle any possible input state. However, deter- 
mining the full process matrix is expensive: it requires in- 
verting at least d 4 — d 2 real parameters for a ^-dimensional 
Hilbert space, in contrast to d 2 parameters for state tomogra- 
phy. Moreover, for analysis of complex molecular dynamics 
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FIG. 1. Analogies between pump-probe and 2D spectroscopy, and 
state and process tomography. By "generalizes," we mean that the 
technique formally contains all information contained in the other. 
By "enables," we mean that the spectroscopic measurement contains 
information that corresponds to a (not necessarily complete) form 
of the specific type of tomography. The top connection, from 2D 
photon-echo to process tomography, was shown by Yuen-Zhou et 
al. 24 . The bottom connection, from pump-probe to state tomography, 
is the focus of this paper. 



in condensed phases, such information, although potentially 
helpful, is not necessary, because most trajectories contained 
in the process matrix start from initial conditions that are im- 
plausible for a molecular aggregate. Indeed, typical theoreti- 
cal investigations of dynamics in light harvesting systems 26 ' 27 
follow dynamics after excitation for only a limited set of plau- 
sible initial states, such as the states which absorb sunlight 
or excitations from neighboring antenna complexes. Finally, 
the relative simplicity of state tomography helps to simplify 
consideration of new theoretical approaches to tomography, 
particularly because process tomography is often based on as 
series of state tomographies 25 . The relationship between these 
prior works and our focus here on pump-probe and state to- 
mography is clarified by the analogies sketched in Figure 1 . 

In this paper, we present a new approach to state tomogra- 
phy of excitonic systems based on pump-probe spectroscopy. 
Our approach is based on a two stage protocol that separates 
the easy (field based) and hard (system based) parts of the 
inversion process. This yields several advantages over prior 
approaches, including the ability to use arbitrarily shaped laser 
pulses and to perform the first inversion even when the second 
inversion is not possible. After presenting the details of each of 
these inversions, we demonstrate their feasibility by applying 
them to invert the simulated spectra of a model dimer with a 
Markovian environment. We close with a consideration of the 
conditions under which state inversion would be feasible for 
a natural light-harvesting system, the FMO complex of green 
sulfur bacteria. 



II. RECIPE FOR PUMP-PROBE SPECTROSCOPY 

We begin by presenting the specific theoretical formal- 
ism for pump-probe spectroscopy that we propose to invert. 
The measured signal in any 3rd-order spectroscopy experi- 
ment, including pump-probe, is a function of the 3rd-order 
polarization 1 . This 3rd-order polarization depends on three 
interactions between the applied fields and the sample, with 



the time-ordering of these interactions enforced by time de- 
lays of the pulses and by looking at the signal emitted in a 
particular phase-matched direction. For a pump-probe exper- 
iment, the first two interactions happen with the same pulse, 
the pump, and the last interaction is with the probe pulse. The 
phase-matched condition is that the signal be observed in the 
direction of the probe. From a formal perspective, this guar- 
antees that the first two interactions (with the pump) happen 
on opposite sides of the density matrix. Based on this phase- 
matched geometry and the response function formalism 13 (see 
Appendix A), we can combine the allowed time orderings to 
write the nonlinear polarization for a pump-probe experiment 
under the rotating wave approximation as 

POO 

P&(t) = / dt 3 Rjv(t 3 ,p$(t-t 3 ))E+(t-t 3 ), (1) 
Jo 

in terms of the pump-probe response 
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This pump-probe response is of identical form to the linear 
response function 1 , but with the electronic ground-state den- 
sity matrix po replaced by the second-order contribution to 

(2) 

the density matrix p pp ; that contributes to the phase-matched 
signal observed in a pump-probe experiment (that is, with 
signal wave- vector fes = fe pr ). The quantities E^{t) and 
E+ (t) denote the complex envelopes of the pump and probe 
pulses, respectively, with E~ = (E + )*. The dipole opera- 
tors = ^2 n d n a n and p^ = (p^)\ with a n as the 
annihilation operator for an electronic excitation on pigment 
n and d n the corresponding dipole moments. The Liouville 
space operator G(t) is the retarded material Green function 
for evolution for time t and • = [p^\ •]. Formally, the 
portion of the second order contribution to the density matrix 
which contributes to the signal is given by 



dt 2 dt 1 G(t 2 )V i±) G(t 1 )V iT) po 



xE±(t-t 2 )E^(t-t 2 -t 1 ). (3) 



In deriving Eqs. (2-3), we employed the rotating wave approx- 
imation and neglected those terms from the double-quantum- 
coherence contribution (fes ^ fepr). Accordingly, we can safely 
neglect the possibility of multiple excitations in the calculation 

(2) 

of p pp ; . In Appendix B we prove that the excited state portion 

(2) 

of pp F J is both equal to the excited state portion of the full 
density matrix and is itself a valid (but unnormalized) density 
matrix. 

The core of our proposed quantum state tomography is the 
sequential inversion of Eqs. (1-3). The remainder of this sec- 
tion discusses additional details relevant to simulating experi- 
mental signals to test our inversion procedure. We emphasize 
that these expressions hold under the very general conditions, 
requiring only the rotating wave approximation, that all ap- 
plied fields are weak and negligible overlap between pump 
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and probe pulses. No assumptions were made concerning the 
shapes of these pump and probe pulses. Our decomposition 
is similar to the window-doorway picture for the pump-probe 
signal 1 , but here we have separated out the influence of the 
control fields, even when not in the impulsive "snapshot" limit. 
Related expressions in terms of a convolution of pump and 
probe components have been shown to facilitate analysis of 
pump-probe experiments with shaped probes 28 . 

There are several obvious extensions or alternatives to this 
recipe that could be useful. For example, the exact same re- 
lations in Eqs. (1-2) hold for a general photon-echo (or non- 
rephasing) experiment with two distinct pump pulses, except 
that in this situation the sum in Eq. (3) to determine the por- 
tion of // 2 ) that contributes to the signal should be removed to 
leave only one of the two phase-matched contributions. Thus 
the photon-echo (PE) signal depends on the second-order den- 
sity matrix given by 

xE+(t-t 2 )E^(t-t 2 -t 1 ). (4) 

The disadvantage of this approach is that these density matri- 
ces do not necessarily contain the total excited state density 
matrix, as is the case for pump-probe. Another option is to 
formally switch the order of the convolution and the response 
function in Eqs. (1-2). In this case, we would be calculating 

(2) 

the pump-probe response of the convolution of p pp ; with the 
probe field. This is easiest to see by considering Eq. (5) be- 
low, and noting that all terms that depend on r, including the 
integral, could be moved to be adjacent to the density matrix 
inside the response function. However, as we shall see, invert- 
ing the convolution is typically the easier step. It therefore 
makes sense to apply the convolution last, so that it can be 
undone first. 



A. Detection scheme and probe convolution 

In a typical pump-probe experiment, the probe pulse has 
a fixed time-envelope, subject to a variable delay time T 
between the two pulses. Accordingly, we may substitute 
E+(t) = E pr (t — T). Likewise, experimental signals are 
most directly interpreted in the frequency domain, so we now 
consider the Fourier transform of the nonlinear polarization, 
P( 3 \uj) = Jdt e iuj ^- T ^P^\t), calculated relative to the 
probe delay T. We can also write the pump-probe response 
in the Fourier domain, i?pp(cj, p^) = J dt e lU)t R?v(t, Pyp). 
In terms of these quantities in the frequency domain with the 
explicit probe delay T, we can then replace Eq. (1) with a 
one-dimensional convolution, 

/CO 
dr R PP (oJ, P$(t))E pi (t - T)e^- T \ 
-co 

(5) 

To obtain this relation, we substituted ts = t — r and extended 
the lower limit of the integral in Eq. (1) to — oo, because by 



definition G(t) = for t < 0. We cannot simply turn this 
convolution into a multiplication by taking the Fourier trans- 
form of these quantities with respect to T, because for small 
or negative delay times T, there are contributions from signals 
where the pump does not necessarily interact before the probe. 
If i^pp(cj, ppp^(r)) does not vary appreciably over the dura- 
tion of the probe pulse, then the equation above is the Fourier 
transform of the probe field envelope, so we can approximate 

P^\u,T) « E pr (uj)R FF (u,pg\T)). (6) 

In the limit of a completely impulsive probe, E pT (t) w EoS(t) 
and thus E pv (uu) is constant, so the nonlinear polarization and 
the pump-probe response are equal up to a constant of propor- 
tionality. 



We measure the non-linear polarization P^ (t) by detect- 
ing the corresponding signal field E$(t) oc iP^ (t) 1 . Here we 
consider heterodyne detection, either with the probe pulse as in 
a standard "self-heterodyned" pump-probe setup, or with a sep- 
arate local-oscillator (LO) pulse. The use of a separate local- 
oscillator is possible in a transient-grating setup, in which the 
pump pulse is replaced by two otherwise identical pumps with 
different wavevectors k\ and k 2 , so that the signal wavevector 
ks = —ki~\-k 2 -\-ks does not match the probe wavevector 
k%. Mathematically, this transient-grating signal yields the 
same non-linear polarization as in pump-probe, although it 
raises experimental complications by requiring phase- stability 
with an additional pulse. In heterodyne detection, the absolute 
value squared of the sum of the signal and local-oscillator (or 
probe) fields can be spectrally dispersed and measured in the 
frequency domain 29 . Typically, the signal field is much smaller 
than than of the local oscillator, so upon subtracting away the 
local oscillator contribution, the measured signal S(lj) is pro- 
portional to Re[Es(u))El (u))], and thus 

S(u>,T) cc lm[P^(u,T)El (uj)]. (7) 

This equation is a multiplication in the frequency domain. 
Hence it is a convolution, and takes on similar form to Eq. (5) 
when expressed in the time-domain. In the pump-probe setup, 
E L0 = E pr , so the signal for a fast probe given by inserting 
Eq. (6) yields 

S(u>,T) a \E pr (uj)\ 2 lmR PF (uj,pg\T)). (8) 

In this case, the signal only depends on the imaginary (absorp- 
tive) part of the non-linear polarization P^ and the pump- 
probe response. In the alternative transient grating setup, as 
long as one is still in the limit of a fast probe, applying a tt/2 
phase shift to the now distinct local oscillator pulse allows for 
obtaining the real (dispersive) part of the pump-probe response 
function in a similarly direct manner 29 . More generally, het- 
erodyne detection with and without a tt/2 phase shift allows 
for obtaining both real and imaginary parts of the non-linear 
polarization, respectively. 
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B. Pump-probe response function 

To isolate the effect of the probe, the pump-probe response 
function given by Eq. (2) can be written as 



i?pp(£, Ppp*) 



Tr 



(9) 



with the pump-probe response operator V(t) defined as 



r(t) 



(-)G(t)V<+> = l -[^-\t),^+\0)}, (10) 



where ^'(t) denotes /i^) in the Heisenberg picture. 
This is similar but not equivalent to a family of quantum 
measurements 4 parametrized by the continuous time variable 
t (or frequency uj in the Fourier domain), since R?? can be 
complex valued. Accordingly, the pump-probe response can 

(2) 

be interpreted as the projection of p pp ; onto V(t), where these 
are viewed as vectors in Liouville space 1 , 



(2) 
PPP 



(ID 



Individual components of the pump-probe response operator 
((V(uj)\a)} are equivalent to the species associated spectra of 
the state \a)) 21 . 

In most spectroscopy experiments, the signal is an ensem- 
ble measurement summed over all possible molecular orienta- 
tions and static disorder of Hamiltonian parameters (inhomo- 
geneous broadening). Accordingly, the pump-probe response 
in Eq. (2) should be replaced by its average over molecular 
orientations and static disorder. The orientational average can 
be handled elegantly using the expression for the pump-probe 
response in Eq. (9): in the magic angle 6 w 54.7° (MA) rel- 
ative polarization configuration between the pump and probe 
pulses 30 , the quantities V(t) and can simply be replaced 
by their independent isotropic averages, 



#fp(*,Pfp) 



= Tr 



MA 



V(t) 



(2) 
PPP 



(12) 



By virtue of the properties of isotropically averaged tensors 31 , 
these independent isotropic averages are equal to the average 
of the quantities obtained from the xx, yy and zz configu- 
rations. In contrast, the ensemble average over static disor- 
der cannot be factorized this way in general, because under 
static disorder the pump-probe operator and density matrix are 
correlated, and altering the system Hamiltonian (e.g., to shift 
transition energies) changes both quantities systematically. 



III. INVERSION PROTOCOLS 

A. Deconvolution of the pump-probe signal 

The first stage of our inversion protocol is a double- 
deconvolution to determine the complex valued pump-probe 

(2) 

response function i? PP (T, p pp ) from the results of a series of 



heterodyne measurements, i.e., the signal S(uj,T). We need 
such a double-deconvolution procedure because the results of 
heterodyne detection depend on a (trivial) convolution over 
the non-linear polarization, which in turn depends on a convo- 
lution over the response function [see Eqs. (5) and (7)]. Since 
the excited state density matrix is entirely contained in the 
pump-probe response function (Appendix B), this inversion 
retains all information about the quantum state. However, it 
is not immediately clear that the real (dispersive) part of the 
response function contains useful information independent of 
the imaginary (absorptive) part, which is the portion measured 
by usual pump-probe experiments. We will return to this point 
in Section V. For now, we make the conservative choice of 
considering inversion of the full complex valued quantity. 

Inverting the signal to obtain the pump-probe response func- 
tion is a non-trivial but important task, since, as pointed out 
above, the signal is directly proportional to the response only 
when the probe pulse is much faster than all energy transfer dy- 
namics. Such pulses can be difficult to realize experimentally. 
The need for a full inversion to obtain the response function 
is particularly relevant for understanding experiments which 
show fast oscillations due to quantum beats, whether these 
are of electronic, vibrational or mixed origin. In such cases, 
the fast probe assumption of Eq. (6) is not valid. We shall 
refer to the use of this approximate description for inversion 
as the "naive" approach. In contrast, a proper treatment of this 
inversion problem should attempt to undo the convolution in 
Eq.(5). 

The convolution of the pump-probe response with the probe 
pulse in Eq. (5) and the non-linear polarization with the local- 
oscillator in Eq. (7) are both particular cases of a Fredholm 
integral equation of the first kind. An extensive literature exists 
on numerical inversion of such equations, known as discrete 
inversion problems 32 . In general, the discretization of such an 
integral equation can be written as 



b = Ax + e, 



(13) 



where b is the measured signal [e.g., the nonlinear polariza- 
tion P 3 (cj,T)], x is the desired quantity to invert [e.g., the 
response function i?pp(cj, T)]. A is a linear operator represent- 
ing the integral equation with appropriate coefficients (deter- 
mined here by the probe field E w ) and e denotes some addi- 
tive experimental error inherent in the data collection. The 
obvious solution to estimating x from A and b is to calculate 
x = A~ x b. However, in practice A may not be invertible. 
Moreover, the presence of even a vanishingly small amount of 
experimental noise e makes exact least- squares minimization 
unsuitable; it will over-fit the noise component e. Accordingly, 
to calculate a robust estimate x of x we use general form 
Tikhonov regularization 32 , 

x = argmin{||Ax - b\\ 2 + A 2 ||Lcc|| 2 } . (14) 

X 

Tikhonov regularization, also known as ridge regression, can 
be derived formally from the perspective of Bayesian in- 
ference, given normally distributed errors and priors 33 . It 
can be equivalently be expressed as the linear least-squares 
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problem, min 1 1 [ £ L ]x — [ o ] 1 1 2 , and thus the exact solution 
is given by ic = [ £ L ] + [ q ] , where + denotes the Penrose- 
Moore pseudoinverse 32 . Ideally, the linear operator L and 
(real) regularization parameter A are chosen so that A 2 1 \Lx\ | 2 
is an optimally weighted penalty on undesirable features of 
the solution x, reflecting our prior knowledge of the gen- 
eral form of x. Common choices of L include the iden- 
tity matrix / and finite-difference approximations to the first 
or second derivative given by (Dix) n = x n — x n _\ and 
(D2x) n = x n+ i — 2x n + x n -\. There are a number of power- 
ful techniques for calculating efficient approximate solutions 
to Eq. (14), especially in cases where the linear operator A is 
structured 34 , such as in our case, where A is a Toeplitz matrix. 

In this paper, we consider a simple, easily parallelized ap- 
proach based on two sequential general form Tikhonov in- 
versions to invert the pump-probe response from heterodyne 
measurements. First, we invert the non- linear polarization 
P( 3 ) (cj, T) from the measured signal S(uj,T) recorded at each 
choice of delay time T. The relationship between these sig- 
nals is simple multiplication by the probe (or local oscillator) 
field in the frequency domain, so this step only uses Tikhonov 
inversion to smooth the reconstruction along the cj-axis. Sec- 

(2) 

ond, we invert the response function iiW(co>, pp P ; (T)) from 
one- dimensional deconvolutions of the non-linear polarization 
P(cj,T), for each fixed value of uj. For this inversion, we 
only use experimental data with the delay between pump and 
probe pulses long enough so that we can ignore pulse overlap 
effects. Otherwise, we would be including non-pump-probe 
contributions to the signal. However, we also reconstruct the 
pump-probe response at shorter times to appropriately handle 
boundary conditions, since the probe convolution means that 
these values for the response function contribute to the nonlin- 
ear polarization inside our region of interest. This results in an 
absence of data for the smallest delay times T, which guaran- 
tees that the convolution matrix A is not invertible, since b has 
a lower dimensionality than x. 

An important additional step is appropriate selection of reg- 
ularization parameters A and L for each of the two inversion 
steps. In practice, there are a wide variety of techniques for 
making these selections and the best choice depends on the 
particular problem at hand 32 . By generating a large number of 
examples of simulated noise, we found the best performance 
for our numerical example of a dimer spectrum (described in 
Sec. IV C) to be obtained with L = D2 (the finite difference 
second derivative operator), and choosing A to minimize the 
generalized cross-validation estimate of the error. The details 
of generalized cross-validation and of our parameter selection 
tests are described in Appendix C. Since in general the re- 
sponse function should have similar types of features to our 
test-example, we expect these choices should work well for 
inverting almost any pump-probe spectra. 

For several hundred time-delays or probe frequencies, we 
find that we can solve each deconvolution with Eq. (14) ex- 
actly and quickly (~1 s on a modern CPU) by calculating the 
singular value decomposition of the matrix [ £ L ] . In principle, 
it would be possible to solve both steps in the inversion of the 
pump-probe response in a single two-dimensional Tikhonov 
regularization. Such 2D inversions are routinely performed in 



image processing 32 , but would require slower, more approxi- 
mate techniques than the exact solution we used here. Since 
we find significant improvement without invoking these more 
complicated methods, we do not use them here. 

This first step in the inversion of pump-probe experiments 
requires only the detection results, i.e., the signal S(lj, T), and 
an excellent characterization of the probe and local oscillator 
fields. No system information is required at all. Likewise, we 
have sacrificed no information from our measurement about 
the internal system information, including its quantum state. 
Thus in principle this stage can be performed with high ac- 
curacy for any system, no matter how complex its internal 
degrees of freedom. 



B. Inverting the quantum state 

The second step to complete the state tomography is to in- 
vert the pump-probe response function Rpp(uj, p^ (T)) to ob- 
(2) 

tain the quantum state pp P ; (T). This is certainly the harder 
step, since it requires the ability to construct the pump-probe 
response of arbitrary states. The necessary information is con- 
tained in the pump-probe operator V{t) given by Eq. (10); cal- 
culating this requires both the transition dipole moments and 
a model for dynamics of the 1-exciton coherences between 
the probe and signal interactions. Essentially the same in- 
formation is necessary to implement proposed algorithms for 
quantum process tomography 23,24 . 

There are two general approaches we could use to obtain 
the necessary system information to perform state tomography: 
(1) the "bottom up" approach of calculating the spectra from an 
independently determined model Hamiltonian, and (2) the "top 
down" approach of estimating the contributions to the pump- 
probe response directly from experimental fits, in the spirit of 
global analysis 21 . These approaches have complementary ad- 
vantages and disadvantages, so in practice some combination 
seems most suitable. For example, the bottom-up approach 
needs careful checking and possible refinement to ensure that 
it consistently describes experimental results. In contrast, the 
top-down approach of global analysis avoids the inherent dif- 
ficulties of modeling dynamics, but needs other input (e.g., 
crystal structures) to relate spectroscopic observables to mi- 
croscopic parameters such as pigment location. In practice, 
neither approach may offer enough information to reliably de- 
code the spectra of complex pigment-protein complexes from 
experimental measurements, since there remain large uncer- 
tainties and discrepancies in the model Hamiltonians of even 
the most extensively studied photo synthetic complexes. We 
will consider these issues in more depth in Section V. 

Here we consider a simple bottom-up protocol for state to- 
mography, based on an assumed model for calculating the 
pump-probe response. It is by no means the only such pos- 
sible state tomography protocol: we choose it because it is 
straightforward to implement, and turns out to be relatively 
robust to imperfections such as static disorder. To perform the 
inversion, we propose to extract an estimate of the excited- state 
electronic density matrix p e (r) from the estimated response 
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function i?p P (u^, r) at that delay time r, for each frequency 
uoi matching the single-exciton transition energies. The rela- 
tionship between the vector of pump-probe response measure- 
ments and the density matrix elements at any fixed time delay 
is linear [see Eq. (11)], so as long as this map is non-singular, 
we can solve for the density matrix by simply applying the ma- 
trix inverse to the vector formed by these estimated response 
function points. It is possible that in some circumstances this 
reconstruction will not yield a valid density matrix, since we 
did not include the constraint that the reconstruction be pos- 
itive semi-definite. In this case, then a best estimate to mini- 
mize the mean-squared-error of the reconstruction should be 
obtained using techniques based on maximum likelihood 35 , al- 
though we do not encounter this issue for the examples we 
consider in this paper. 

An additional important practical step is the choice of Liou- 
ville space in which the extracted state lies. Our results so far 
hold for transition dipole operators and time evolution without 
any particular restrictions concerning electronic vs vibrational 
states or the Hamiltonian we use to describe our system. How- 
ever, as a practical matter for excitonic energy transfer in light 
harvesting systems, we are most interested in inverting the 

(2) 

electronic degree of freedom. The electronic portion of pp P ; 
has useful structure: namely, it only includes nonzero elements 
in two blocks, the 0- and 1 -excitation subspaces. We denote 
the projection of the density matrix p onto these subspaces by 
p g and p e , respectively. There is only one electronic state in 
the 0-excitation subspace (the ground state g), so the electronic 
portion of p g must be in that state, \g) (g\. In the Markov limit, 
or for delay times much longer than the bath relaxation time, 
the vibrational portion of p g will be in thermal equilibrium 
Pg q . These facts determine p g , up to a constant of proportion- 
ality: the ground state population. Because total probability 

(2) 

is conserved in the process of laser excitation, Tr pp P ; = 0, 
so the ground state population is related to the excited state 



population by Tr pg 



(2) 



(2) 

Tr p K e J . Accordingly, we can write 



o( 2 ) 



(15) 



In this case, the pump-dependence in the pump-probe signal 
[Eq. (3)] is entirely contained in the excited state portion of 

(2) 

the density matrix. Since for weak fields p e w p e , with 
Eqs. (11) and (15) we have a linear map from any excited state 
density matrix p e to the corresponding pump-probe response. 
For a system with n electronic states, we can parameterize this 
unnormalized density matrix in terms of a linear combination 
of n 2 real parameters 5 , since the excited state density matrix 
is positive (see Appendix B) and thus Hermitian. 

A brief discussion of the scalability of this approach is in 
order. In this work we focus on inverting the electronic por- 
tion of the excited state of a molecular aggregate in an open 
environment with exactly one relevant electronic excitation 
per pigment molecule. Based on the real and imaginary parts 
of the pump-probe response function in the magic angle con- 
figuration, our state tomography protocol in principle has 2n 
independent real parameters with which to extract the n 2 real 
parameters (including normalization) necessary to describe an 



arbitrary excited state density matrix of n pigments 5 . Accord- 
ingly, we cannot necessarily expect this procedure to scale 
beyond a dimer (n = 2), for which we numerically demon- 
strate the success and stability of this inversion procedure in 
the next section. The recently proposed quantum process to- 
mography algorithm based on peak and cross-peak amplitudes 
in 2D spectroscopy 24 has similar scaling difficulties. In gen- 
eral, process tomography requires determining n 4 — n 2 real 
parameters in the process matrix from at most 12n 2 possible 
experimental measurements: the real and imaginary signals, n 
coherence and n rephasing frequencies, at most 3 independent 
polarization configurations and 2 phase-matched geometries. 
These estimates, however, hold only for this specific approach 
and with a randomly oriented ensemble. For example, there 
are 9 different potential probe/signal- field polarization config- 
urations for an oriented sample, so with such measurements 
we could potentially invert systems as large as n = 18. Indeed, 
in Section V, we discuss how generalizations of this procedure 
could nonetheless be used to successfully invert spectra for sin- 
gle complexes of the Fenna-Matthews- Olson complex (FMO), 
with seven molecular pigments. 



IV. EXAMPLE: DIMER MODEL 

To understand in more detail how the quantum state de- 
termines the pump-probe signal, we consider the case of the 
signal for a dimer of coupled pigments. In general, we can 
write an effective Hamiltonian for the electronic excited states 
of a dimer in the form 

H Q \ — E\a\ai + E 2 a\a2 + J{a\ai + a[a2). (16) 

The terms E\ and E 2 are the transition energies of sites 1 and 2, 
and J is the pigment-pigment coupling energy. We restrict the 
system to at most one excitation on each site, so our state space 
is spanned by the set {\g), |ei), |e2), |/)}, denoting the ground 
state, excitation of the first or second site, and excitation of 
both sites. We further assume the usual linear coupling to 
a bath of phonons. Details of the bath are specified below. 
The electronic part of this Hamiltonian can be diagonalized by 
applying a unitary rotation U to the single- excitation subspace, 
given by 



U 



cos sin 
— sin 6 cos 6 



(17) 



where we defined the mixing angle 6 = \ arctan(2J/A) with 
A = Ei — E 2 . These single excitation eigenstates are denoted 
\a) and The transition dipole moments for each pigment 
are di and d 2 , oriented with relative angle (j). 



A. Analytical calculation of pump-probe response 

To begin, we choose a parametrization for the excited state 
density matrix of a dimer. In general, any valid density matrix 
for a two-level system can be written in any basis in form 
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\ (I + r • a), in terms of the Pauli matrices tr = {cr x , o- y ,(J z } 
and the Bloch vector r = {ri, r2, ^3}, with real and |r| 2 < 
l 4 . This can be straightforwardly generalized to unnormalized 
density matrices by adding the normalization ro and defining 
<To = L in which case the set of valid but unnormalized states 
are those that can be written in the form \ (r • cr), where r is 
now the 4-dimensional vector {ro, ri, r2, r%} with r\ + r| + 
^3 < ^0 an d r o > 0- We will use these four real parameters 
to parametrize the excited state electronic density matrix p e 
for our tomography protocol, since it has population ro C 1. 
Using this representation, the total second-order correction to 
the electronic density matrix from Eq. (15) is given by 



P y w J = -ro\g)(g\ + ■ 



(18) 



with ^ = diCti, where di is the component of the dipole- 
transition vector parallel to the probe polarization. For conve- 
nience, we define f a and f' a to denote the Fourier transform 
of the time evolution operator that leads to a peak in the pump- 
probe spectrum at frequency u a , with decay constant 7 or 7', 
where the prime indicates the decay constant for the transition 
between the 1- and 2-exciton manifolds instead of between the 
0- and 1-exciton manifolds: 



fa 



1 



i(u a - cu) - 7' 



f 

J a 



i(u a - lj) - 7' 



(20) 



We define fp and fa analogously, for the components peaked 
at uop. 



For convenience, we suppose the state is written in terms of 
the eigenbasis expansion of the excited states { \a) , so the 
parameters r\ and r2 correspond to excitonic coherences and 
rs corresponds to the balance of population between excitonic 
states. In practice, the experimental signal is only known up 
to a constant factor, so we can only hope to be able to reliably 
determine the normalized excited-state density matrix, given 
by the usual Bloch- vector elements {ri /ro , r2 /ro , r3 /ro } . 

It is now straightforward (if tedious) to write down the ex- 
act pump-probe response function in terms of microscopic 
parameters. For illustrative purposes, we do so here for a 
dimer with a Markovian bath described by Redfield theory 
in the secular approximation 3 . The time-evolution contained 
directly in the pump-probe response function is for times fol- 
lowing the probe interaction, so the relevant part of the sys- 
tem density matrix for this evolution includes coherences be- 
tween ground and singly excited states and between singly 
and doubly excited states. In secular-Redfield theory, coher- 
ences in the excitonic basis only evolve with exponential decay, 
G(T)\a)(b\ = e _7abT 0(T)|a)(6|, where G(T) is the retarded 
material Green function denoting evolution for time T, G is 
the Heaviside step function, and 7 a 6 is some complex number 
with positive real part. Since the formulas for the response 
function will accordingly be most compact in the exciton ba- 
sis, we consider the excitonic transition dipole moments given 

by, 



l^gcx = Mi cos + P2 sin 
p g p = —pi sin 6 + P2 cos 
l^af — Mi sin + P2 cos 
Pf3f = p\ cos — P2 sin (9, 



(19a) 
(19b) 
(19c) 
(19d) 



Using Eq. (11), the calculation of the pump-probe response 
function for an arbitrary state is determined by the vector- 
ized version of the pump-probe operator, ((V(u)\. For this 
dimer problem, we define the pump-probe bra vector such that 

Rff = ((P\r)}, where |r)) = r = {r , ri, r 2 , r 3 }. Such a 

(2) 

relation still holds upon substitution of |r)) for pp P \ since the 
relation between the two given by Eq. (18) is linear. With this 
convention, evaluating the pump-probe operator in Eq. (10) for 
this dimer model described by secular-Redfield theory yields 
the general result, 



«7>|oc 



(fa + fp) " 
Vga^g(3 (fa fp) 

Wpff'ot 



^gafa 



(fa + f'p) 

(fa-f'p) 
~ V 2 g pfp + Vaff'p 



(21) 



This equation holds for each single molecule that would con- 
tribute to the pump-probe signal. We can also calculate the 
exact isotropic average of Eq. (21) over an ensemble of ran- 
domly oriented molecules. In terms of the original Hamilto- 
nian parameters, it is given by 



(cos 2 6 + 5 2 sin 2 0) (/; - 3/„) + (sin 2 + 5 2 cos 2 0) (f' p - 3f p ) + 5 sin20 C os^ (-f' a +f' p - 3/ Q + 3/„)l 
-\ {5 2 - 1) sin20 (f' a + f'p + f a + fp) + <5cos20cos0 [f' a +/£-/*- fp) 
i [\ (S 2 - 1) sin20 (f' a -fp + f a - fy) + ,5cos20cos0 (-f' a + /£ + /„- /^)] 
- (cos 2 + 5 2 sin 2 0) (/; + f a ) + (sin 2 + S 2 cos 2 6) (f' p + f p ) + S sin 20 cos <f> (f' a + f' f3 -f a - fp) 

(22) 
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where 6 is the excitonic mixing angle, S = \ d<i \ / \ d\ | is the 
ratio of the two site transition dipole moments and <p is the 
angle between them. Note that neither of these equations in- 
cludes the effects of static disorder, which could be accounted 
for by averaging the pump-probe response function over each 
member of the ensemble. Formally, it does not suffice to sep- 
arately average ((P|, since under static disorder the state |r)) 
also varies in correlated way (see Sec. II B). 

Equation (22) makes it possible to determine some cases in 
which inverting the isotropically averaged state cannot possi- 
bly be successful, regardless of the exact inversion protocol. 
We can identify these cases because successful inversion re- 
quires that the elements of (((V\)i so be linearly independent. 
For example, in the homodimer case with both pigments fixed 
to have the same transition energies (6 = 7r/4 or 6 = 3tt/4) 
and equal transition-dipole moment magnitudes (a = 1), the 
pump-probe signal does not depend on the coherence terms, 
so it will be impossible to invert them (n and r2). Likewise, 
the coherence terms do not contribute if the transition dipole 
moments have identical magnitude (S = 1), and either are 
oriented perpendicularly (cos (j) = 0) or there are matching 
dephasing rates for the 0-1 and 1-2 coherences (f a = f a 
and fp = fL as occurs in the high-temperature limit). As 
Yuen-Zhou et al. found for the same dimer model 23 , quantum 
process tomography also fails under similar but not identical 
conditions. 



B. Numerical example 

For a numerical example, we consider the dimer model 
used in a prior investigation of quantum process tomography 23 . 
We model excitation by a resonant 40 fs full-width-at-half- 
maximum (FWHM) pump centered at 12 800 cm -1 . The pa- 
rameters in the electronic Hamiltonian are E\ — 12 881 cm -1 , 
E2 = 12 719cm -1 and J = 120 cm -1 , and the experiment 
is performed on an ensemble with normally distributed static 
disorder of standard deviation 40 cm -1 added to each site en- 
ergy. The transition dipole moments are fixed with ratio S = 
|^2/^i| = 2 and orientation angle <p = 0.3. Each pigment 
is assumed to be coupled to an independent bath of phonons, 
with spectral density of the form J(uj) = ^-uje - ^^ with 
uo c = 120 cm -1 and A = 30 cm -1 . The bath is assumed to 
be at thermal equilibrium at T = 273 K and is modeled by 
secular Redfield theory 3 , including only the real (dissipative) 
part of the Redfield tensor. 

To simulate an experimental dataset, we first calculate the 
non-linear polarization P^(cj,T) for a probe of the same 
shape as the pump pulse, on a grid of 181 probe frequen- 
cies uj (intervals of 3.33 cm -1 between 12 500 cm -1 and 
13 100 cm -1 ) and 140 central time-delays T between pump 
and probe pulses (intervals of 6.81 fs between 50 fs and 1 ps). 
From this non-linear polarization, we then calculate the results 
of a hypothetical heterodyne detection with a local oscilla- 
tor matching the probe pulse, with and without a tt/2 phase 
shift. Finally, we accounted for noise in detection by including 
additive noise with uniformly random phase and amplitude 
drawn from a standard deviation with width equal to 10 -2 
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FIG. 2. Absorptive (left) and dispersive (right) parts of the pump- 
probe response function i?pp(cj, (r)) (top) and the corresponding 
heterodyne detected signal S(uj,t) (bottom) for our dimer model 
system. The dashed line indicates the two exciton transition energies 
in this system. Only the absorptive part (left) is revealed directly 
by a pump-probe experiment. Obtaining the dispersive part (right) 
requires a transient grating setup with heterodyne detection, as de- 
scribed in Sec. II A. 



(high noise) or 10 - 6 (low noise) times the maximum ampli- 
tude over all delay times and frequencies of the heterodyne 
detected signal S(u,T). These simulated measurements (for 
high noise) are shown in Fig 2, together the response function 
from which they are calculated. Note that we chose this low 
noise threshold of 10 -6 because we found that our response 
function inversion algorithm ran into numerical errors when 
the inversion was performed with less added noise, because 
the optimal damping for the Tikhonov regularization in these 
cases is at A = 0. Fortunately, this should not be a serious 
concern for actual experimental measurements, for which it is 
always possible to add a small amount of simulated noise. 



C. Response function inversion 

Figure 3 illustrates the performance of the completed dou- 
ble Tikhonov regularization based deconvolution algorithm for 
typical low and high noise examples of our test problem. We 
compare with the "naive" approach of assuming that the probe 
is impulsive and using Eq. (6) to obtain the response function 
by simply dividing the signal by absolute value squared of the 
probe field \E pr (uj)\ 2 . Table I summarizes the results of 1000 
such simulated inversions for each of the high and low noise 
examples. In addition to the double Tikhonov and naive meth- 
ods, we also consider the alternatives of substituting the naive 
approach individually for each of the two Tikhonov steps. Re- 
call that in the first stage of the inversion (S — » P^), the 
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FIG. 3. (a) Example reconstruction of the pump-probe response at 
fixed probe-frequency uj a for an instance of the high-noise test prob- 
lem, (b) Errors in the inverted pump-probe response obtained by the 
direct and Tikhonov inversion methods for a single example of the 
low and high noise test problems. The error is given by the abso- 
lute value squared of the difference between the estimated and actual 
response function, \R pp (uj, r) — i?pp(cj, r)| 2 . 
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TABLE I. Summary of deconvolution performance over 1000 in- 
stances of simulated experimental noise. Improvement is the mul- 
tiple of the reduction in mean-squared-error compared to the direct 
approach. Uncertainties, indicating one standard deviation in the em- 
pirical distribution, are unmarked for the 10 ~ 6 noise level since there 
is essentially no observed variation in the inversion between different 
noise instances. 



Tikhonov regularization serves only to smooth the data. In the 
second stage (P^ — ^ i?pp), the Tikhonov regularization also 
performs a deconvolution over the probe envelope. 

The results in Fig. 3 and Table I show that the Tikhonov 
based inversion is a clear improvement over the naive ap- 
proach. The low and high noise examples allow us to observe 
that the Tikhonov regularizations remove two types of errors 
inherent in the naive inversion: (1) errors from noisy measure- 
ments and (2) errors associated with the convolution of the 
pump-probe response over the finite probe duration. In the 
high noise case, both errors are large; in the low noise case, 
the second source of error dominates. Clearly, reducing the 
experimental noise associated with measurement alone does 
not suffice to accurately estimate the pump probe response: 
for example, the double Tikhonov inversion under high noise 
outperforms the naive inversion under low noise by a factor of 
7 in the mean-squared-error (MSE), despite a reduction in the 
level of noise by four orders of magnitude. 

As Figure 3 shows, the errors in the estimate of the response 
function are not uniformly distributed, revealing structure rele- 
vant to our specific example and also to more generic systems. 
The largest errors are associated with smallest delay times. 
This makes sense, since the smallest delay times are those at 
which at the response function (shown in Fig. 2) varies most 
rapidly. For almost any system, the pump-probe spectrum 
will change fastest at short delay times, but this is especially 
true for our example system, where the pump-probe signal in- 
cludes contributions from quickly oscillating coherences. The 
Tikhonov estimates face an additional stability challenge at 
short delays times, since, as discussed above, the reconstruc- 
tion cannot use measurements from the pulse overlap regime. 



D. State tomography 

Since we have demonstrated that the first, response func- 
tion inversion can be performed with vanishing error, we now 
consider inverting the exact pump-probe response function to 
obtain the state of our model dimer. Despite the presence of 
static disorder in our example, we use the factorization of the 
response function in Eq. (12). We are obliged to do so even 
though strictly speaking the relationship does not hold, be- 
cause the alternative of reconstructing the density matrix for 
each member of the ensemble from a bulk measurement is un- 
realistic. Accordingly, even without adding noise associated 
with the measurement, when carried out for an ensemble, our 
inversion faces potential stability issues because of the static 
disorder. However, for a single complex, this inversion would 
be robust. 

The degeneracies and near degeneracies in Eq. (22) mean 
that for most Hamiltonian choices our inversion algorithm can 
only robustly extract at most three of the four parameters nec- 
essary to fully characterize the dimer excited state, since the 
reconstruction matrix will be poorly conditioned. The condi- 
tion number of a linear transformation gives a bound on the 
multiplicative increase in the relative error after performing 
the linear transformation 34 . For our specific numerical exam- 
ple, the condition number drops from 3700 to 3.1 when we 
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include only three parameters. One source of these stability is- 
sues for a dimer is evident from Eq. (22): since our numerical 
examples has well separated transition energies, the main con- 
tribution to the peaks in the dispersive part of the signal is to 
the imaginary part of the coherence term, T2 . This leaves our 
inversion to recover three parameters (7*0, T\ and rs) from the 
two peak amplitudes in the absorptive signal. Since recovering 
three unknowns from two equations is not possible, we need 
to fix one of these values in order to make the inversion stable. 
The obvious choice is to fix the normalization 7*0, since the 
total excited state population should remain constant after the 
end of the pump pulse until spontaneous decay, on timescales 
approaching 1 ns for natural pigment-protein complexes 2 . To 
determine the normalization, we solve for it at a moderately 
long delay time (e.g., r = 10 ps) at which point we can safely 
assume (at least under secular Redfield dynamics) that the real 
part of the coherence 7*1 —> 0, but very few excitations have 
been lost. If these timescales are not easily separable, then this 
normalization term could be fit to an exponential decay. 

The results of applying this state tomography procedure to 
our numerical example with varying levels of static disorder 
are shown in Fig. 4. The fidelity, ranging from to 1, provides 
a numerical summary of the quality of the state reconstruction 4 . 
For the level of static disorder chosen by Yuen-Zhou et al. 
(40 cm -1 ), the reconstruction [panel (a)] has a worst-case fi- 
delity of 99.5% over delay times T shorter than 1 ps, and an 
average-case fidelity of 99.9%. However, we see that both 
the worst-case and average-case fidelities drop sharply as the 
static disorder is increased above this level [panel (b)], since 
our assumption that the pump-probe response can be factor- 
ized between the pump-probe projection and the second order 
density matrix becomes increasingly unrealistic. 



V. SCALING TO LARGER SYSTEMS 

Can state tomography scale to systems larger than an ex- 
citonic dimer? In particular, can we apply it to precisely re- 
veal the excitonic state in a natural light-harvesting system? 
Any scaling difficulties will be encountered in the second step 
of our inversion protocol, to invert the quantum state from 
the response function, since the relationship between the re- 
sponse function and measured signal does not directly depend 
on system parameters. Successful state tomography certainly 
requires both knowledge of how each density matrix element 
contributes to measurements and appropriate conditions such 
that each element, at least in principle, makes an independent 
and non-zero contribution. By construction, these conditions 
were satisfied for our hypothetical dimer example. We found 
that the primary limitation on inverting the state was ensemble 
disorder, which can in principle be avoided by single molecule 
techniques. Now, to explore the limits of state tomography, we 
relax these assumptions in order to consider the feasibility of 
state tomography in an actual protein-pigment complex. 

As a model light-harvesting system, we focus on a monomer 
of the Fenna-Matthews-Olson (FMO) complex of green sul- 
fur bacteria, which consists of 7 pigment molecules 18 ' 36 . The 
FMO complex is a widely used model system for understand- 
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FIG. 4. Results of quantum state tomography for our dimer test prob- 
lem, (a) Original (solid) and reconstructed (dotted) values for each 
element of the Bloch state vector for the reconstruction with static 
disorder of standard deviation 40 cm -1 . Normalization is omitted 
since the state vector elements are rescaled such that ro = 1 fixed 
for all times following initial excitation, (b) Worst- and average-case 
fidelities for the reconstructions p e (r) for delay times r in the range 
50 fs to 1 ps as a function of the width (standard deviation) of the 
distribution of static disorder. Results are obtained from an ensemble 
average over 10 6 samples for each point. 



ing photosynthetic energy transfer and is thus is one of the 
best characterized natural protein-pigment complex. The crys- 
tal structure for the FMO complex is known, which combined 
with input from spectroscopy experiments, has allowed for 
general agreement on an electronic Hamiltonian 36,37 . Because 
the arrangement of pigments is fairly disordered, each excited 
state in a monomer of the FMO complex is bright, although 
they overlap in the presence of homogeneous and inhomoge- 
neous broadening. This is important, since full state tomogra- 
phy would certainly not be possible on a system with multiple 
dark states, because no optical probes could reveal the distri- 
bution of energy among these states. However, typical of the 
situation for other natural pigment-protein systems, there is 
little consensus on the magnitude of the static disorder or the 
spectral density of the electronic-vibrational coupling 27 . These 
difficulties are compounded by the theoretical and computa- 
tional challenge of modeling dynamics in a system as large as 
FMO exactly for arbitrary system-bath interaction strength 26 . 
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For our concrete example, we use the electronic Hamiltonian 
for FMO of Chlorobaculum tepidum from Ref. 36, with the 
spectral density and computational model of secular Redfield 
theory matching those used in for the dimer example. Note 
that this model includes Gaussian distributed static disorder 
with standard deviation 42.5 cm" 1 (100 cm" 1 FWHM). 

Our formalism for the pump-probe response function allows 
us to place bounds on the feasibility of any state tomography 
procedure, since the relationship between system information 
and the resulting pump-probe spectra is entirely contained in 
the pump-probe response operator V(u). By looking at the 
ensemble average of this operator, we implicitly consider in- 
version under the scenario that the average over static disorder 
can be factorized between the pump-probe operator and the 
density matrix. This assumption was shown to be success- 
ful in the dimer example above when the magnitude of static 
disorder was not too large. As discussed in Section II B, the 
pump-probe operator at each frequency can be interpreted as a 
Liouville space bra- vector ((V(u)\. Accordingly, it is possible 
to interpret the calculation of a pump-probe response as the 
act of applying the linear operator 

P = J duj\uj)((V(oj)\ (23) 

(2) 

to the state |pp P ; )) . We now consider the properties of the linear 
operator P in the limit of effectively continuous sampling of 
probe frequencies uj. To represent states in Liouville space, 
we use a basis set that allows us to represent each state with 
n 2 real values, in terms of populations \n)(n\ and coherences 
\n)(m\ + \m)(n\ and i\n)(m\ — i\m)(n\. This allows us to 
construct a real- valued version of the map P that takes real 
valued state vectors to real valued spectra by concatenating the 
real and imaginary parts of P. 

To begin, in Figure 5 we plot the elements of the absorp- 
tive (real) part of the pump-probe operator for our model of 
the FMO complex in the isotropic average (magic angle con- 
figuration). We represent the operator in terms of the species 
associated spectra of each exciton population and the real and 
imaginary part of each coherence, so that the pump-probe spec- 
tra of any excited state is equal to the linear combination of 
the plotted spectra weighted by the indicated density matrix 
elements. In addition to the unperturbed spectra, we also plot 
the range of possible spectra given current uncertainty about 
the best fit parameters. We conservatively estimate the uncer- 
tainty in the electronic Hamiltonian by sampling over additive 
independent Gaussian noise of width 20 cm" 1 for each site 
energy and 10% of the value of each off-diagonal coupling. 
This uncertainty is in addition to the static disorder, which we 
leave with fixed magnitude. The most striking feature of these 
spectra is that, at least in the isotropic average, the dominant 
contribution to the pump-probe spectra is from the population 
terms. The smaller contribution of most coherence terms, com- 
pounded by the already smaller values of the coherences in the 
density matrix due to dephasing, explains why it is difficult 
to observe oscillations due to electronic coherence in pump- 
probe spectra 18 . Even for extremely precise measurements, 
the uncertainty in some of these species associated spectra 



suggests that our current Hamiltonian characterization does 
not suffice to reliably obtain most density matrix elements. 
Indeed, the dominance of the diagonal terms suggests that a 
practical scheme for partial state tomography could consist of 
entirely ignoring the off-diagonal terms. 

Another approach to estimating the feasibility of inversion 
for arbitrary states is to look at the spectral properties of the 
operator P as revealed by the singular value decomposition, 
P = USV\ where U and V are unitary and S is diagonal 
with positive elements. In particular, we focus on the singular 
values di, given by the diagonal elements of S in descend- 
ing order and normalized to the highest singular value g\ . To 
compare the feasibility of inversion under various conditions, 
we plot these singular values in Figure 6. The singular values 
reveal significant information about the feasibility of an inver- 
sion: in general, inversion is more feasible when the singular 
values Gi decay more slowly 32 . For example, the condition 
number, which gives an upper bound on the ratio by which the 
relative error can increase in an inversion, is equal to the ratio 
of the largest to the smallest singular values. In Figure 6, the 
condition number is one over the value shown for i = 49. 

Figure 6 clearly indicates the relative significance of differ- 
ent experimental constraints insofar as they affect state tomog- 
raphy. We see two major changes that reduce the condition 
number of the inversion by around two orders of magnitude 
each. First, reducing temperature from 300 K to 77 K improves 
the conditioning because spectral features in the species asso- 
ciated spectra are sharpened, as shown in Figure 6. Second, 
changing from a randomly oriented (isotropic) to an oriented 
sample helps because in principle 9 times more independent 
measurements are possible than in the magic angle configura- 
tion, one for each combination of x, y and z polarizations for 
the probe and signal interactions. This is illustrated in Figure 7, 
which plots the peak amplitudes of each of the species associ- 
ated spectra. In some cross-polarization configurations, coher- 
ences give contributions of similar magnitude to populations. 
The plot of singular values also suggests that the dispersive 
component of the signal, which as discussed in Section II A 
requires a transient grating experiment, offers little additional 
information, at least with this particular model, compared to 
the absorptive component, which can be obtained from pump 
probe directly. However, an independent local-oscillator pulse 
would still be necessary to separately resolve cross polariza- 
tions between the probe and signal fields from an oriented 
sample. 

We now consider the improvement in performance that is 
offered by oriented samples and by single molecule (i.e., 'sin- 
gle complex') experiments. While current experimental state 
of the art is limited to measurement of non-oriented samples, 
several technological developments suggest that pump-probe 
experiments on oriented light harvesting systems may be pos- 
sible in the foreseeable future. For example, it was recently 
demonstrated that it is possible to perform ultrafast spectro- 
scopies on crystallized proteins 38 . The possibilities offered by 
single molecule spectroscopic techniques are even more sig- 
nificant for improving the performance of the inversion. An 
important limitation of the above analysis based on singular 
values is that it does not explicitly examine the failure of the 
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FIG. 5. Species associated spectra, defined by the contribution of the marked density matrix elements to the pump-probe response, for the FMO 
complex at 77 K (blue) and 300 K (red), obtained as the average of 1000 samplings over static disorder. Labels indicate the contributing density 
matrix element in the excitonic basis. Shaded regions indicate central 95% confidence intervals obtained from 1000 additional samplings over 
Hamiltonian uncertainty, as described in the text. 



factorization between pump-probe operator and density oper- 
ator in the case of an ensemble measurement with static dis- 
order. As Figure 4(b) showed for the dimer problem, this can 
be a serious issue if there is significant variation between dif- 
ferent members of the ensemble. Resolving such difficulties 
inherent in ensemble measurements would become possible 
with experiments on single molecules. For single molecule 



measurements, which in some sense are also inherently on ori- 
ented samples, there is no constructive interference across the 
sample to allow for satisfying a phase-matching condition, so 
it is not possible to directly measure the non-linear polariza- 
tion. However, essentially equivalent information can be ob- 
tained from a measurement of the non-linear florescence with 
phase-cycling 39 . Although the statistics for single molecule 
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FIG. 6. Normalized singular values from a singular value decompo- 
sition of the real valued pump-probe map P given by Eq. (23) for 
the FMO complex under various conditions. The top line is from the 
spectra of a single monomer at 77 K, from the combination of mea- 
surements in all independent polarization configuration. Subsequent 
lines add additional constraints, which apply cumulatively: ensemble 
measurement (over static disorder), the isotropic average of the signal, 
only the absorptive (real) part of the signal and finally performing the 
measurement at room temperature. 
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FIG. 7. Maximum amplitude over probe frequencies of the species 
associated spectra for an FMO monomer at 77 K for (a) the isotropic 
average and (b) each independent polarization configuration of the 
probe and local oscillator, including the ensemble average over static 
disorder. The labeling of each species matches that used in Figure 5: 
all entries including and above the diagonal correspond to the real 
part of the matching density matrix element (in the excitonic basis), 
and all entries below the diagonal correspond to the imaginary part. 
The cartesian coordinates were chosen arbitrarily, matching those 
used in an assignment of the crystal structure. 



measurements would assuredly be worse, these measurements 
do offer the best conditioned inversions, as indicated by the 
plot of singular values in Fig. 6. 



VI. CONCLUSIONS 

In this paper we have presented a new approach to the in- 
verse problem of extracting the excited state of a molecular 
aggregate via ultrafast spectroscopy. The theoretical analysis 
of this approach and its demonstrations made here for a dimer 
and for FMO suggest a number of insights about how to de- 
sign ultrafast spectroscopy experiments and the limits of the 
information we can obtain from them. 

For example, a particular benefit of our approach is that 
it is well suited to evaluating the benefits of employing col- 
ored pulses or an ultrafast pulse shaper, because our formal- 
ism made no assumptions about the shape of pump and probe 
pulses: they are only required not to overlap in time. In con- 
trast, prior proposals for process tomography relied on the 
assumption that interactions with laser pulses are much faster 
than the timescale of dissipative system dynamics 23,24 . Ac- 
cordingly, we envision potentially using our scheme for verifi- 
cation of state preparation following shaped pump pulses 40 ' 41 , 
scenarios for which single molecules or oriented samples are 
similarly helpful. Indeed, pump-probe spectroscopy has been 
used to experimentally verify ultrafast coherent control 42 . The 
deconvolution step in our inversion protocol also made no as- 
sumptions about the shape of the probe pulse. Although our 
extensive numerical simulations found no cases in which a 
shaped probe pulse was preferable to the corresponding time- 
frequency limited pulse, our inversion protocol can just as 
easily invert the non-linear response from, say, pump-probe 
measurements where the probe has residual chirp. This sug- 
gests the possibility of improving pulse characterization for 
better time resolution instead of or in addition to efforts to 
further compress pulses in time 28 . Likewise, it is straightfor- 
ward to invert the results of a series of measurements with 
different shaped or colored probe pulses: unless it is needed 
for time-resolution, a broadband probe could be replaced by 
a set of different colored probes that together span all probe 
frequencies of interest. 

Our multi-stage inversion also has clear extensions to more 
general non-linear spectroscopies beyond QST and pump- 
probe. As we discussed in Sec. II, we could apply essen- 
tially the same inversion procedure to a photon-echo exper- 
iment (or in the non-rephasing geometry) to determine the 
phase-matched component of the 2nd-order density matrix that 
contributes to the observed signal [Eq. (4)]. If it is possible 
to construct a full set of independent phase-matched initial 
conditions, our state tomography procedure could be used for 
process tomography, along the lines of a previous proposal 23 . 
More generally, the same type of two stage inversion suggests 
a new indirect approach for process tomography: first, field 
information should be used to invert the 3rd-order response 
function; second, system information should be used to extract 
the process matrix. This approach of inverting the response 
function as an intermediate quantity allows us to be sure we 
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have obtained the maximum information from experiments be- 
fore considering the harder theoretical problem of extracting 
system parameters, such as the state, process matrix or under- 
lying Hamiltonian. Such extensions will be pursued in future 
work. 

Finally, we point out that we demonstrated our inversion 
approach in this work mostly for ensemble systems, with av- 
eraging over molecular orientations and static disorder. These 
features inevitably reduce the performance of the inversion. 
However, with the ongoing advances in realizing spectroscopy 
on oriented samples and on singe molecule systems, it seems 
not unrealistic to expect high fidelity inversion of pump-probe 
and other 2D spectroscopies in the near future. 
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Appendix A: Derivation of pump-probe signal 

In third-order spectroscopy, the phase-matched signal de- 
pends on the third-order polarization 13 , which can be written 

as 



JJJ dt 3 dt 2 dh R { k}{hMM) 

xE^(t-t 3 )E^(t-t 3 -t 2 ) 
xE^(t-t 3 -t 2 -h) 



(Al) 



assuming that the interactions happen in the numbered or- 
der and invoking the rotating wave approximation, which 
is excellent for these systems. The quantity k$ = uiki + 
u 2 k 2 + ^3&3 is the signal wave- vector, (ui,u 2 ,u 3 ) £ 
{(— , +, +), (+, — , +), (+, +, — )} correspond to the three ex- 
perimental geometries with non-zero signal (rephasing, non- 
rephasing and double-quantum-coherence, respectively) and 
E+(t) denotes the complex profile of the ith pulse (we use 
the convention E~ = (£? + )*). The system dynamics are con- 
tained in the phase-matched components of the third order 
response function (£3, t 2 , ti), which is given by 



Rkht 3 ,h,h) - 



Tr 



^G^V^G^V^G^V^po 



(A2) 



in terms of the quantities defined in Sec. II. 

In a pump-probe experiment, the first two interactions are 
with the same pulse (k\ = k 2 ), called the pump, and the signal 
is observed in the direction ks = £3, so u\ = — u 2 . The third 
interaction is with the probe field. Accordingly, the signal 
is given by adding together the rephasing and non-rephasing 



interactions, and both pulse orderings 1-2-3 and 2-1-3. Re- 
arranging the terms that result from inserting Eq. (A2) into 
Eq. (Al) yields Eqs. (1-3). 



Appendix B: Positivity of /4 2) 

We consider the system density matrix p in the presence of 
weak fields of strength e«l. Let p^ denote the contribu- 
tion to the density matrix of strength 0(e n ). Thus we can write 
p = p (0) + + p (2) + p (3) + 0(e 4 ). Since the temperature 
is typically several orders of magnitude smaller than the elec- 
tronic energy gap, the system starts in a tensor product of the 
electronic ground state \g)(g\ and the equilibrium vibrational 
stote^p™ = \g)(g\®f% { . 

(2) 

We are interested in the excited state portion of pp P \ the 
component of p^ that contributes to the phase-matched sig- 
nal ks = fcpr given by Eq. (3). We write this excited state 

(2) (2) 

density portion as p\' = Qpp P , where Q denotes the pro- 
jection onto the 1 -excitation manifold. Since the non-phase- 
matched components involve 2-excitation states, we also have 
pi^ = Qp^ '. This projection is given by Qp = I\pl\, where 
^1 = Em \ m )( m \ is me identity operator restricted to the 
1 -excitation manifold and \m) denotes the state where only 
pigment m is excited. Since the map Q is written in the appro- 
priate form and E n ^i^n = ^> 2 is completely positive, with 
< Tr Qp < l 4 . 

Moving from the ground state to the 1 -excitation manifold 
requires at least two applications of the creation/annihilation 
operators contained in the dipole transition operators and 
a dipole operator must be applied an even number of times. 
Accordingly, Qp^ = for n = 1, 2,4, which leaves Qp = 
pi^ + 0(e 4 ). Since we proved Qp is positive and p^ itself 
is 0(e 2 ), we have shown that p^ is positive, up to relative 
errors of 0(e 2 ). 



Appendix C: Parameter selection for Tikhonov 
regularization 

To perform deconvolution using Eq. (14), we need practi- 
cal procedures for determining L and A. Since there are no 
universal best choices for all inverse problems 32 , we took the 
approach of testing the performance of different techniques on 
simulations matching the dimer problem we analyze in Sec- 
tion IV. To begin, we compared the performance of general 
form Tikhonov regularization with L equal to /, D\ and D 2 , 
with A chosen optimally so as to minimize the exact mean- 

squared-error \\x — x\\ 2 . A summary of reconstructions of 

(2) 

Rpp(uj, pp P ; (T)) for 1000 instances of low and high noise is 
shown in Table II. 

As a benchmark, we consider the ratio of the mean- squared- 
error from the Tikhonov estimate to the mean- square-error of 
our naive estimate R(uj, T) ex P(cj, T)/E w (oo), which holds 
in the limit of an instantaneous probe pulse [Eq. (6)]. For 
our state inversion algorithm, reconstruction of the response 
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Noise 


Penalty 


Selection 


^opt 


Improvement 
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Exact 


0.151 ±0.007 


2.3 ±0.2 


10" 


-2 




Exact 


0.340 ± 0.034 


4.9 ±0.7 


10" 


-2 


D 2 


Exact 


0.67 ±0.11 


6.2 ±1.2 


10" 


-2 


D 2 


GCV 


0.78 ±0.17 


5.7 ±1.2 


10" 


-2 


D 2 


NCP 


2.57 ±0.36 


2.3 ±0.4 


10" 


-3 


I 


Exact 


0.049 ±0.017 


5.9 ±0.5 


10" 


-3 


D x 


Exact 


0.073 ±0.010 


30 ±5 


10" 


-3 


D 2 


Exact 


0.112 ±0.017 


56 ± 12 


10" 


-3 


D 2 


GCV 


0.103 ±0.019 


52 ±12 


10" 


-3 


D 2 


NCP 


0.449 ± 0.039 


19 ±2 



TABLE II. Regularization performance for different penalty operators 
and parameter selection techniques for 1000 instances of random 
noise with relative magnitude 10 _2 orl0 -3 . Numbers are the mean 
plus or minus one standard deviation. Improvement is the multiple of 
the reduction in mean-squared-error for the reconstructed response 
function using Tikhonov regularization over the error associated with 
the naive impulse-probe estimate. 



function is most important at frequencies matching the exci- 
ton transition energies, so we picked u = uj a , the transition 
frequency of the higher energy exciton state. We found qualita- 
tively similar results for uj = ujp and other choices of uj as well. 
As Table II shows, with exact selection of the best regulariza- 
tion parameter A, we found consistently best performance with 
D 2 , the linear operator approximating the second derivative of 
the response function with respect to the delay time T. This 
is an intuitively reasonable choice, since plausible response 
functions should be smooth. 

With the choice for L determined, we also need a realistic 
procedure for selecting the regularization parameter A. In a 
true inversion problem the response function x is unknown, 
so we cannot choose A to minimize the exact mean-squared- 
error. There are a variety of standard techniques for making 
this choice, with performance that can vary widely depending 
on the problem being solved, so selection of an appropriate 
method requires more tests on simulated data. We considered 
two such methods for L = D 2 : generalized cross-validation 
(GCV) and the normalized cumulative periodogram (NCP). 
We calculated the GCV error using the exact pseudoinverse 
solution 33 and the NCP error by adding together the errors for 
the real and imaginary parts of the spectra 32 . The results, also 
shown in Table II, show that GCV is the best choice for our 
test problem, with performance nearly matching that of exact 
selection technique. In contrast, NCP systematically overes- 
timated the noise, as indicated by regularization parameters 
about four times larger than the exact selection method. How- 
ever, NCP still offered an improvement in the reduction of the 
mean-squared-error compared to the naive approach. 
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